﻿#pragma once

#include <learn/geom/geom_basis.h>

namespace xi
{

namespace geom
{

template <std::size_t N, bool Rat = false> class Geom_BSplineSurface : public Geom_BasisSurface<N>
{
    static constexpr std::size_t D = N + Rat;
    using VecR = VecT<D, double>; // (wP,w)
    using Vec = VecT<N, double>;  // P

  public:
    static constexpr bool Rational = Rat;
    Geom_BSplineSurface(Geom_BSplineSurface &&) = delete;

    /**
     * @brief 基本构造函数
     *
     * @param[in] cpts
     * @param[in] uKnots
     * @param[in] vKnots
     * @param[in] uDeg
     * @param[in] vDeg
     * @param[in] normalize
     */
    Geom_BSplineSurface(std::vector<VecR> &&cpts, std::vector<double> &&uKnots, std::vector<double> &&vKnots,
                        std::size_t uDeg, std::size_t vDeg, bool normalize = true)
        : m_uDeg(uDeg), m_vDeg(vDeg), m_uKnots(std::move(uKnots)), m_vKnots(std::move(vKnots)), m_cpts(std::move(cpts))
    {
        ASSERT_INFO(uDeg == u_degree() && vDeg == v_degree(), "Degree not match.");
        ASSERT_INFO(m_cpts.size() == u_pole_num() * v_pole_num(), "Control points not match.");

        // 如果为有理 Bezier 曲面，则对点的坐标进行规范化 (P,w) -> (wP,w)
        if constexpr (Rat)
        {
            if (normalize)
            {
                for (auto &p : m_cpts)
                {
                    double w = p[N];
                    p[N] = 1.0;
                    p *= w;
                }
            }
        }
    }

    /**
     * @brief u 方向的次数
     *
     * @return std::size_t
     */
    std::size_t u_degree() const
    {
        return m_uDeg;
    }

    /**
     * @brief v 方向的次数
     *
     * @return std::size_t
     */
    std::size_t v_degree() const
    {
        return m_vDeg;
    }

    /**
     * @brief u 方向的 knot 数
     *
     * @return std::size_t
     */
    std::size_t u_knot_num() const
    {
        return m_uKnots.size();
    }

    /**
     * @brief v 方向的 knot 数
     *
     * @return std::size_t
     */
    std::size_t v_knot_num() const
    {
        return m_vKnots.size();
    }

    /**
     * @brief 获得 u 方向的 knot
     *
     * @param[in] i
     * @return double
     */
    double u_knot(std::size_t i) const
    {
        return m_uKnots[i];
    }

    /**
     * @brief 获得 v 方向的 knot
     *
     * @param[in] i
     * @return double
     */
    double v_knot(std::size_t i) const
    {
        return m_vKnots[i];
    }

    /**
     * @brief 设置 u 方向的 knot
     *
     * @param[in] i
     * @param[in] k
     */
    void set_u_knot(std::size_t i, double k)
    {
        m_uKnots[i] = k;
    }

    /**
     * @brief 设置 v 方向的 knot
     *
     * @param[in] i
     * @param[in] k
     */
    void set_v_knot(std::size_t i, double k)
    {
        m_vKnots[i] = k;
    }

    /**
     * @brief u 控制点个数
     *
     * @return std::size_t
     */
    std::size_t u_pole_num() const override
    {
        return u_knot_num() - m_uDeg - 1;
    }

    /**
     * @brief v 控制点个数
     *
     * @return std::size_t
     */
    std::size_t v_pole_num() const override
    {
        return v_knot_num() - m_vDeg - 1;
    }

    /**
     * @brief 获得控制点
     *
     * @param[in] u
     * @param[in] v
     * @return Vec
     */
    Vec pole(std::size_t u, std::size_t v) const override
    {
        if constexpr (Rat)
            return homogenous_to_cartesian(m_cpts[u * v_pole_num() + v]);
        else
            return m_cpts[u * v_pole_num() + v];
    }

    /**
     * @brief 设置控制点
     *
     * @param[in] u
     * @param[in] v
     * @param[in] p
     */
    void set_pole(std::size_t u, std::size_t v, const Vec &p) override
    {
        if constexpr (Rat)
            m_cpts[u * v_pole_num() + v] = cartesian_to_homogenous(p) * weight(u, v);
        else
            m_cpts[u * v_pole_num() + v] = p;
    }

    /**
     * @brief 获得权重
     *
     * @param[in] u
     * @param[in] v
     * @return double
     */
    double weight(std::size_t u, std::size_t v) const
    {
        if constexpr (Rat)
            return m_cpts[u * v_pole_num() + v][N];
        else
            return 1.0;
    }

    /**
     * @brief 设置权重
     *
     * @param[in] u
     * @param[in] v
     * @param[in] w
     */
    void set_weight(std::size_t u, std::size_t v, double w)
    {
        static_assert(Rat, "Rational surface not supported.");
        m_cpts[u * v_pole_num() + v] *= w / weight(u, v);
    }

    /**
     * @brief 获得无穷远点
     *
     * @param[in] u
     * @param[in] v
     * @return VecR
     */
    VecR inf_pole(std::size_t u, std::size_t v) const
    {
        static_assert(Rat, "Rational surface required.");
        return m_cpts[u * v_pole_num() + v];
    }

    /**
     * @brief 设置无穷远点
     *
     * @param[in] u
     * @param[in] v
     * @param[in] p
     */
    void set_inf_pole(std::size_t u, std::size_t v, const VecR &p)
    {
        static_assert(Rat, "Rational surface required.");
        m_cpts[u * v_pole_num() + v] = p;
    }

    /**
     * @brief 获得 u 节点
     *
     * @return const std::vector<double>&
     */
    const std::vector<double> &u_knots() const
    {
        return m_uKnots;
    }

    /**
     * @brief 获得 v 节点
     *
     * @return const std::vector<double>&
     */
    const std::vector<double> &v_knots() const
    {
        return m_vKnots;
    }

    /**
     * @brief 移动 u 节点
     *
     * @return std::vector<double>
     */
    std::vector<double> move_u_knots()
    {
        return std::move(m_uKnots);
    }

    /**
     * @brief 移动 v 节点
     *
     * @return std::vector<double>
     */
    std::vector<double> move_v_knots()
    {
        return std::move(m_vKnots);
    }

    /**
     * @brief 获得控制点
     *
     * @return const std::vector<VecR>&
     */
    const std::vector<VecR> &poles() const
    {
        return m_cpts;
    }

    /**
     * @brief 移动控制点
     *
     * @return std::vector<VecR>
     */
    std::vector<VecR> move_poles()
    {
        return std::move(m_cpts);
    }

    /**
     * @brief 参数范围
     *
     * @return std::array<double, 4>
     */
    std::array<double, 4> bounds() const override
    {
        return {m_uKnots[u_degree()], m_uKnots[m_uKnots.size() - u_degree() - 1], m_vKnots[v_degree()],
                m_vKnots[m_vKnots.size() - v_degree() - 1]};
    }

    /**
     * @brief 计算曲面上点的位置
     *
     * @param[in] u
     * @param[in] v
     * @return Vec
     */
    Vec d0(double u, double v) const override
    {
        auto p = surface_point(m_cpts.data(), m_uKnots.data(), m_vKnots.data(), m_uKnots.size(), m_vKnots.size(),
                               u_degree(), v_degree(), u, v);
        if constexpr (Rat)
            return homogenous_to_cartesian(p);
        else
            return p;
    }

    /**
     * @brief 计算曲面上点的切向量
     *
     * @param[in] u
     * @param[in] v
     * @return std::array<Vec, 3>
     */
    std::array<Vec, 3> d1(double u, double v) const override
    {
        if constexpr (Rat)
        {
            auto ders = rat_surface_deriv(m_cpts.data(), m_uKnots.data(), m_vKnots.data(), m_uKnots.size(),
                                          m_vKnots.size(), u_degree(), v_degree(), 1, u, v);
            return {ders[0][0], ders[1][0], ders[0][1]};
        }
        else
        {
            auto ders = surface_deriv_alg1(m_cpts.data(), m_uKnots.data(), m_vKnots.data(), m_uKnots.size(),
                                           m_vKnots.size(), u_degree(), v_degree(), 1, u, v);
            return {ders[0][0], ders[1][0], ders[0][1]};
        }
    }

    /**
     * @brief 计算曲面上点的二阶导数
     *
     * @param[in] u
     * @param[in] v
     * @return std::array<Vec, 6>
     */
    std::array<Vec, 6> d2(double u, double v) const override
    {
        if constexpr (Rat)
        {
            auto ders = rat_surface_deriv(m_cpts.data(), m_uKnots.data(), m_vKnots.data(), m_uKnots.size(),
                                          m_vKnots.size(), u_degree(), v_degree(), 2, u, v);
            return {ders[0][0], ders[1][0], ders[0][1], ders[2][0], ders[0][2], ders[1][1]};
        }
        else
        {
            auto ders = surface_deriv_alg1(m_cpts.data(), m_uKnots.data(), m_vKnots.data(), m_uKnots.size(),
                                           m_vKnots.size(), u_degree(), v_degree(), 2, u, v);
            return {ders[0][0], ders[1][0], ders[0][1], ders[2][0], ders[0][2], ders[1][1]};
        }
    }

    /**
     * @brief 序列化
     *
     * @return utils::Json
     */
    utils::Json dump_json() const override
    {
        utils::Json json;
        json["type"] = "bpsline_surface";
        json["rational"] = Rat;
        json["uDegree"] = int(u_degree());
        json["vDegree"] = int(v_degree());
        json["uKnots"] = m_uKnots;
        json["vKnots"] = m_vKnots;
        json["poles"] = m_cpts;
        return json;
    }

    /**
     * @brief u 方向插入 r 次
     *
     * @param[in] u
     * @param[in] r
     * @return true
     * @return false
     */
    bool insert_u_knot(double u, std::size_t r = 1)
    {
        auto span = find_span(m_uKnots.data(), m_uKnots.size(), u_degree(), u);
        auto s = find_mult_span(m_uKnots.data(), span, u_degree(), u);
        if (s + r > u_degree())
            return false;

        std::vector<VecR> Q((u_pole_num() + r) * v_pole_num());
        std::vector<double> U(u_knot_num() + r);
        surface_knot_ins<VecR, 0>(m_cpts.data(), m_uKnots.data(), u_knot_num(), v_knot_num(), u_degree(), v_degree(), u,
                                  span, s, r, Q.data(), U.data());

        m_uKnots = std::move(U);
        m_cpts = std::move(Q);

        return true;
    }

    /**
     * @brief v 方向插入 r 次
     *
     * @param[in] v
     * @param[in] r
     * @return true
     * @return false
     */
    bool insert_v_knot(double v, std::size_t r = 1)
    {
        auto span = find_span(m_vKnots.data(), m_vKnots.size(), v_degree(), v);
        auto s = find_mult_span(m_vKnots.data(), span, v_degree(), v);
        if (s + r > v_degree())
            return false;

        std::vector<VecR> Q(u_pole_num() * (v_pole_num() + r));
        std::vector<double> V(v_knot_num() + r);
        surface_knot_ins<VecR, 1>(m_cpts.data(), m_vKnots.data(), u_knot_num(), v_knot_num(), u_degree(), v_degree(), v,
                                  span, s, r, Q.data(), V.data());

        m_vKnots = std::move(V);
        m_cpts = std::move(Q);

        return true;
    }

    /**
     * @brief 细化 u 方向
     *
     * @param[in] x
     * @return true
     * @return false
     */
    bool refine_u_knot(std::vector<double> &&x)
    {
        if (x.empty())
            return false;

        std::vector<VecR> Q((u_pole_num() + x.size()) * v_pole_num());
        std::vector<double> U(u_knot_num() + x.size());

        surface_knot_refine<VecR, 0>(m_cpts.data(), m_uKnots.data(), u_knot_num(), v_knot_num(), u_degree(), v_degree(),
                                     x.data(), x.size(), Q.data(), U.data());

        m_uKnots = std::move(U);
        m_cpts = std::move(Q);

        return true;
    }

    /**
     * @brief 细化 v 方向
     *
     * @param[in] x
     * @return true
     * @return false
     */
    bool refine_v_knot(std::vector<double> &&x)
    {
        if (x.empty())
            return false;

        std::vector<VecR> Q(u_pole_num() * (v_pole_num() + x.size()));
        std::vector<double> V(v_knot_num() + x.size());

        surface_knot_refine<VecR, 1>(m_cpts.data(), m_vKnots.data(), u_knot_num(), v_knot_num(), u_degree(), v_degree(),
                                     x.data(), x.size(), Q.data(), V.data());

        m_vKnots = std::move(V);
        m_cpts = std::move(Q);

        return true;
    }

    /**
     * @brief 在 u,v 方向上升次到 n,m 次
     *
     * @param[in] n
     * @param[in] m
     * @return std::tuple<bool, bool>
     */
    std::tuple<bool, bool> elevate(std::size_t n, std::size_t m)
    {
        std::tuple<bool, bool> ret = {false, false};
        n = std::max(n, u_degree());
        m = std::max(m, v_degree());

        std::vector<double> U, V;
        std::vector<VecR> Q;

        // v 方向升阶
        std::size_t row = v_knot_num();
        if (m > v_degree())
        {
            V.resize(v_knot_num() * m);
            Q.resize(V.size() - m - 1);
            row = curve_elevate(m_cpts.data(), m_vKnots.data(), m_vKnots.size(), v_degree(), m - v_degree(), Q.data(),
                                V.data());
            V.resize(row);

            row -= m + 1;
            Q.resize(row * u_pole_num());
            for (std::size_t j = 0; j < u_pole_num(); ++j)
                curve_elevate(m_cpts.data() + j * v_pole_num(), m_vKnots.data(), m_vKnots.size(), v_degree(),
                              m - v_degree(), Q.data() + j * row, V.data());

            m_vDeg = m;
            m_vKnots = std::move(V);
            m_cpts = std::move(Q);
            std::get<1>(ret) = true;
        }

        // u 方向升阶
        std::size_t col = u_knot_num();
        if (n > u_degree())
        {
            U.resize(u_knot_num() * n);
            Q.resize(U.size() - n - 1);
            col = curve_elevate(m_cpts.data(), m_uKnots.data(), m_uKnots.size(), u_degree(), n - u_degree(), Q.data(),
                                U.data(), v_pole_num());
            U.resize(col);

            col -= n + 1;
            Q.resize(v_pole_num() * col);
            for (std::size_t i = 0; i < v_pole_num(); ++i)
                curve_elevate(m_cpts.data() + i, m_uKnots.data(), m_uKnots.size(), u_degree(), n - u_degree(),
                              Q.data() + i, U.data(), v_pole_num(), v_pole_num());

            m_uDeg = n;
            m_uKnots = std::move(U);
            m_cpts = std::move(Q);
            std::get<0>(ret) = true;
        }

        return ret;
    }

  protected:
    std::size_t m_uDeg;
    std::size_t m_vDeg;
    std::vector<double> m_uKnots;
    std::vector<double> m_vKnots;
    std::vector<VecR> m_cpts;
};

} // namespace geom

} // namespace xi
